import sys
import math
import numpy as np
from PyQt5.Qt import *
from numpy import Inf
from numpy.fft import fft,ifft,fft2,ifft2,ifftshift,fftshift
from matplotlib.backends.backend_template import FigureCanvas
from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
from mpl_toolkits.mplot3d import Axes3D
from PyQt5.uic import loadUi
import matplotlib.pyplot as plt
from PyQt5.QtWidgets import *


class Form(QWidget):
    def __init__(self):
        super(Form, self).__init__()
        loadUi('seasurface.ui', self)
        self.figure = plt.figure()
        self.canvas = FigureCanvas(self.figure)
        self.axes = self.figure.add_subplot(111, projection='3d')
        self.verticalLayout_9.addWidget(self.canvas)
#----------JS----------↓----------------------------------------
    @pyqtSlot()
    def on_pushButtonDRAWSEA_clicked(self):
    # def seasurfaceJS(self):
        pi=math.pi
        sin=np.sin
        cos=np.cos
        linspace=np.linspace

        nx = int(self.lineEditSEAX.text().split(',')[2]); xmin =float(self.lineEditSEAX.text().split(',')[0]); xmax = float(self.lineEditSEAX.text().split(',')[1]); x = linspace(xmin,xmax,nx)
        ny = int(self.lineEditSEAY.text().split(',')[2]); ymin =float(self.lineEditSEAY.text().split(',')[0]); ymax = float(self.lineEditSEAY.text().split(',')[1]); y = linspace(ymin,ymax,ny)
        wind_dataU =float(self.lineEditSEAU.text()); wind_datathetaU = 0; wind_dataX = float(self.lineEditSEAUX.text())

        U = float(wind_dataU)
        thetaU = float(wind_datathetaU)

        B = np.fix((U / 0.836)**(2 / 3))
        imunit = 1j

        g = 9.80665
        gxg = g * g
        UxU = U * U
        alpha = 8.1e-3

        dx = x[1] - x[0]
        kxmax = 1 / (2 * dx)
        kx = linspace(-kxmax, kxmax, nx)
        dy = y[1]- y[0]
        kymax = 1 / (2 * dy)
        ky = linspace(-kymax, kymax, ny)
        Kx, Ky = np.meshgrid(kx, ky)
        K = (Kx **2 + Ky **2)**0.5
        F = ((g * K)**0.5) / (2 * pi)
        F[F == 0] = Inf
        K[K == 0] = Inf
        dFdK = ((g/ K)**0.5) / (4 * pi)
        OMEGA = 2 * pi * F

        X = float(wind_dataX) # % 风驱长度或风程
        OMEGAp = 7 * pi * (g/U) * (g*X / U*U)**(-0.33)
        # % 谱峰频率
        fm = OMEGAp / (2 * pi)
        Tp = 2 * pi / OMEGAp
        cgamma = 3.3
        # % 谱峰提升因子，均值为3.3
        csigma0 = 0.07
        if OMEGA.any() > OMEGAp:
            csigma0= 0.09
        # % 当频率（OMEGA）> OMEGAp, csigma0要取0.09，反之取0.07
        calpha = 0.076 * (g * X / U)**(-0.22)
        # % 无因次常数，通常即alfa
        delta = np.exp(-((OMEGA - OMEGAp)** 2)/ (2 * csigma0**2 * OMEGAp**2))
        # % 指数项
        SF = calpha * g*g * OMEGA**(-5)* (np.exp(-(5 / 4) * (OMEGA / OMEGAp)**(-4)))* (cgamma**(delta))
        SK = SF* dFdK

        #扩散谱，。待写。。。。。。。

        white_noise = np.random.randint(-127, 127, (ny, nx))/127
        WHITE_NOISE = fft2(white_noise)
        NOISE_amplitude = abs(WHITE_NOISE)

        NOISE_energy = sum(WHITE_NOISE**2 )
        WHITE_NOISE = WHITE_NOISE / NOISE_energy

        centered_WHITE_NOISE = fftshift(WHITE_NOISE)
        NOISE_phase = np.angle(centered_WHITE_NOISE)

        NOISE_amplitude = NOISE_amplitude* SK

        # % 随机变量的引入
        NOISE_ipart = NOISE_amplitude * sin(NOISE_phase)
        NOISE_rpart = NOISE_amplitude * cos(NOISE_phase)

        filtered_NOISE = NOISE_rpart + imunit * NOISE_ipart
        filtered_NOISE = fftshift(filtered_NOISE)

        # % 基于逆傅里叶变换获取2D粗粗糙面
        s = ifft2(filtered_NOISE)
        s = s.real
        XX,YY=np.meshgrid(x,y)

        plt.clf()
        self.axes.clear()
        ax = Axes3D(self.figure,auto_add_to_figure=False)
        self.figure.add_axes(ax)
        ax.plot_surface(XX, YY, s)
        self.canvas.draw()
#JS-----------------↑------------------------------------

#-----机翼————--------------------↓-----------
    @pyqtSlot()
    def on_pushButtonjiyi_clicked(self):
        plt.clf()
        def camber_line(x, m, p, c):
            return np.where((x >= 0) & (x <= (c * p)),
                            m * (x / np.power(p, 2)) * (2.0 * p - (x / c)),
                            m * ((c - x) / np.power(1 - p, 2)) * (1.0 + (x / c) - 2.0 * p))
        def dyc_over_dx(x, m, p, c):
            return np.where((x >= 0) & (x <= (c * p)),
                            ((2.0 * m) / np.power(p, 2)) * (p - x / c),
                            ((2.0 * m) / np.power(1 - p, 2)) * (p - x / c))
        def thickness(x, t, c):
            term1 = 0.2969 * (np.sqrt(x / c))
            term2 = -0.1260 * (x / c)
            term3 = -0.3516 * np.power(x / c, 2)
            term4 = 0.2843 * np.power(x / c, 3)
            term5 = -0.1015 * np.power(x / c, 4)
            return 5 * t * c * (term1 + term2 + term3 + term4 + term5)
        def naca4(x, m, p, t, c=1):
            dyc_dx = dyc_over_dx(x, m, p, c)
            th = np.arctan(dyc_dx)
            yt = thickness(x, t, c)
            yc = camber_line(x, m, p, c)
            return ((x - yt * np.sin(th), yc + yt * np.cos(th)),
                    (x + yt * np.sin(th), yc - yt * np.cos(th)))
        # naca2412
        # m = 0.02
        # p = 0.4
        # t = 0.12
        # c = 1.0
        m = float(self.lineEditm.text())
        p = float(self.lineEditp.text())
        t = float(self.lineEditt.text())
        c = float(self.lineEditc.text())
        x = np.linspace(0, 1, 200)
        for item in naca4(x, m, p, t, c):
            plt.plot(item[0], item[1], 'b')
            plt.axis('equal')
            plt.xlim((-0.05, 1.05))
            self.axes.plot(item[0], item[1], 'b')
            self.canvas.draw()


if __name__ == "__main__":
    app = QApplication(sys.argv)
    form = Form()
    form.show()
    sys.exit(app.exec_())


